HNOBac Manuscript
  1. Methods
  2. LDH
  • Introduction
  • Methods
    • RNASeq
    • Cell Counts
    • CFUs
    • LDH
    • Cytokines
    • HEK-Blue
  • References
  • R Session Info

Table of contents

  • Data Input and Selection
    • File Paths
    • Data clean-up
    • Saving files
  • Stats and Plots
    • File Paths
    • Function for each temp condition and endpoint
    • Apply to each temp

LDH

library(tidyverse)
library(readxl)
library(ggpubr)
library(ggtext)
library(scales)
library(lme4)
library(afex)
library(emmeans)
library(gridExtra)

Data Input and Selection

File Paths

# Folder paths
input_path <- "data/input_data/LDH/"
metadata_path <- "data/metadata/LDH"

# Create subfolders for output files
dataframes_folder <- "data/dataframes"
if (!file.exists("data/dataframes")) {
  dir.create("data/dataframes", recursive = TRUE)
}

# Load data and metadata
input_data <- read_excel(file.path(input_path, "LAK24_0118_LDHData.xlsx")) %>% mutate_if(is.character, factor)
Bac_order <- read_csv(file.path(metadata_path, "Order_Bacteria_LDH.csv"))

Data clean-up

# Filter only samples matched to CFU wells
LDH_data <- input_data %>% 
  filter(CFUmatched == "T")

# Factor Ordering and Styling
LDH_data <- merge(LDH_data, Bac_order, by = "Species") 
LDH_data$bacteria_label <- factor(LDH_data$bacteria_label, levels = Bac_order$bacteria_label)
LDH_data$Line <- fct_recode(LDH_data$Line, "HNO918" = "A", "HNO204" = "B", "HNO919" = "C") 

# Average technical replicates for each individual experiment
LDH_avg <- LDH_data %>%
  group_by(Date, Line, Species, bacteria_label, Temp, Well_Endpoint, Collection_Time) %>% 
  summarise(avg_Value = mean(Value)) %>%
  # Add labels to distinguish final collection times from time 0
  mutate(Collection_Label = if_else(Collection_Time == 24 | Collection_Time == 48, "final", "initial"))

# Calculate log2 fold change between final time point vs. time 0
LDH_FC <- LDH_avg %>%
  group_by(Date, Line, Species, bacteria_label, Temp, Well_Endpoint) %>%
  reframe(FC = avg_Value[Collection_Label == "final"]/avg_Value[Collection_Label == "initial"]) 

Saving files

# Save data frames as CSV files in the output folder
write_csv(LDH_data, file.path(dataframes_folder, "LDH_values.csv"))
write_csv(LDH_FC, file.path(dataframes_folder, "LDH_FC.csv"))

# Save data frames as R objects in the output folder
saveRDS(LDH_data, file.path(dataframes_folder, "LDH_values.rds"))
saveRDS(LDH_FC, file.path(dataframes_folder, "LDH_FC.rds"))

# Cleaning-up all objects from the environment
rm(list = ls())

# Use this to read the final objects
LDH_data <- readRDS("data/dataframes/LDH_values.rds")
LDH_FC <- readRDS("data/dataframes/LDH_FC.rds")

Stats and Plots

File Paths

# Folder paths
dataframes_path <- "data/dataframes"
metadata_path <- "data/metadata/LDH"

# Create subfolders for output files
figures_folder <- "data/outputs/LDH/figures"
if (!file.exists("data/outputs/LDH/figures")) {
  dir.create("data/outputs/LDH/figures", recursive = TRUE)
}
stats_folder <- "data/outputs/LDH/stats"
if (!file.exists("data/outputs/LDH/stats")) {
  dir.create("data/outputs/LDH/stats", recursive = TRUE)
}

# Load data and metadata
LDH_FC <- readRDS("data/dataframes/LDH_FC.rds")
Bac_order <- read_csv(file.path(metadata_path, "Order_Bacteria_LDH.csv"))

Function for each temp condition and endpoint

# Function to analyze each temp condition
analysis_function <- function(data, each_temp, each_endpoint, cutoff_pvalue, cutoff_FC) {
  
  # Subset the data to the selected temp
  data_subset <- data %>%
    filter(Temp == each_temp) %>%
    filter(Well_Endpoint == each_endpoint)
  
  # Mixed-effects model with random effects
  model <- lmer(FC ~ Species 
                + (1|Line) + (1|Line:Date), 
                data = data_subset)
  #Anova
  anova <- anova(model)
  anova_df <- as.data.frame(anova) %>%
    mutate(sign = case_when(
      `Pr(>F)` < cutoff_pvalue ~ "*",
      TRUE ~ "")) %>%
    mutate_if(is.numeric, ~ format(., digits = 2, scientific = TRUE))
  
  # Calculate Individual contrasts
  emmeans_model <- emmeans(model, ~ Species)
  emmeans_contrasts <- pairs(emmeans_model, adjust = "none")    
  
  # Extract random effects and convert to dataframe (if not singular)
  random_effects_df <- as.data.frame(VarCorr(model)) %>%
    mutate(proportion = round(100 * (vcov / sum(vcov)), 2)) 
  
  # Adds predictions based on fixed effects, averaged over random effects. It gives a population estimate
  data_subset <- cbind(data_subset, predval = predict(model,re.form = NA, se.fit = TRUE))
  data_summary_df <- data_subset %>%
    group_by(Species, bacteria_label) %>%
    summarize(mean.real = mean(FC),
              mean.predval = mean(predval.fit), 
              mean.predval.se = mean(predval.se.fit)) %>%
    mutate(max = mean.predval + 2*mean.predval.se,
           min = mean.predval - 2*mean.predval.se)
  
  # Convert contrasts to dataframe and adjust pvalues. Filter contrast to Uncolonized only
  contrasts_df <- as.data.frame(summary(emmeans_contrasts)) %>%
    filter(str_detect(contrast, "Uncolonized")) %>%
    mutate(p.adj.holm = p.adjust(p.value, method = "holm")) %>%
    mutate(sign = case_when(
      p.adj.holm < 0.05 ~ "*",
      TRUE ~ ""))
  
  # Edits to the contrast dataframe to include pvalue brackets in plot
  contrasts_df <- contrasts_df %>%
    separate(contrast, into = c("condition1", "condition2"), sep = " - ", remove = F) 
  
  contrasts_df$group1 <- Bac_order$bacteria_label[match(contrasts_df$condition1, Bac_order$Species)]
  contrasts_df$group2 <- Bac_order$bacteria_label[match(contrasts_df$condition2, Bac_order$Species)]
  
  # Calculate fold-change values for each contrast
  contrasts_df <- contrasts_df %>%
    ungroup() %>%
    left_join(select(data_summary_df, Species, mean.predval), by = join_by(condition1 == Species)) %>%
    left_join(select(data_summary_df, Species, mean.predval), by = join_by(condition2 == Species), suffix = c(".1", ".2")) %>%
    mutate(FC = mean.predval.1 / mean.predval.2,
           FC = if_else(FC < 1, -1 / FC, FC),
           highlighted = case_when(
             FC <= -cutoff_FC ~ "-",
             FC >= cutoff_FC ~ "+",
             TRUE ~ "")) 
  
  # Select p values to plot and define their location
  contrast_sign <- contrasts_df %>%
    filter(sign != "" & highlighted != "") %>%
    mutate(p.adj.holm = format(p.adj.holm, digits = 2, scientific = TRUE))
  
  location <- max(data_subset$FC, na.rm = TRUE) * 1.1
  
  # Standard Boxplot
  plot_1 <- ggplot() +
    geom_boxplot(data = data_subset, aes(x = bacteria_label, y = FC, fill = bacteria_label)) + 
    
    geom_jitter(data = data_subset, aes(x = bacteria_label, y = FC, shape = Line), 
                fill = "grey50", color = "grey30", size = 2, width = 0.05, stroke = 0.75) +
    
    scale_fill_manual(values = c("grey80","#AA35E3","#2e67f2","#927ed1")) +
    scale_shape_manual(values = c(21, 22, 24)) +
    
    labs(x = "",
         y = "Fold change in LDH from -1 hour") +
    
    theme_bw() +
    theme(panel.grid = element_blank(),
          legend.position = "none",
          text = element_text(size = 20), 
          axis.text.x = element_markdown(), 
          axis.text.y = element_text(color = "black"))
  
  # Plot with predicted means and standard errors of the estimates
  plot_2 <- ggplot() +
    geom_point(data = data_subset, 
                aes(x = bacteria_label, y = FC, fill = bacteria_label, color = bacteria_label, shape = Line), 
                position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
                size = 1.5, alpha = 0.75, stroke = 0.75) +
    
    geom_point(data = data_summary_df, aes(x = bacteria_label, y = mean.predval), shape = 3, size = 3) +
    geom_errorbar(data = data_summary_df, aes(x = bacteria_label,
                                              y = mean.predval,
                                              ymax = max,
                                              ymin = min),
                  width = 0.4) +
    
    scale_fill_manual(values = c("#5b5b5b","#AA35E3","#2e67f2","#927ed1")) +
    scale_color_manual(values = c("#5b5b5b","#AA35E3","#2e67f2","#927ed1")) +
    scale_shape_manual(values = c(21, 22, 24)) +
    
    scale_y_continuous(expand = c(0.1,0)) +
    
    labs(x = "",
         y = "Fold change in LDH from -1 hour",
         fill = "Bacteria", color = "Bacteria", shape = "HNO Line") +
    
    theme_bw() +
    theme(panel.grid = element_blank(), 
          legend.text = element_markdown(),
          text = element_text(size = 20), 
          axis.text.x = element_markdown(), 
          axis.text.y = element_text(color = "black"))
  
  # Conditionally add p-value annotations layer
  if (nrow(contrast_sign) > 0) {
    plot_2 <- plot_2 +
      stat_pvalue_manual(contrast_sign, label = "p.adj.holm", y.position = location,
                         tip.length = 0.02, bracket.shorten = 0.2, vjust = -0.2, bracket.size = 0.3, size = 3.5)
  } else {
    plot_2 <- plot_2
  }
  
  # Arrange plot and tables for summary pdf
  table <- contrasts_df %>%
    select(condition1, condition2, p.adj.holm, sign, mean.predval.1, mean.predval.2, FC, highlighted) %>%
    mutate(p.adj.holm = format(p.adj.holm, digits = 2, scientific = TRUE))
  
  Tmin <- ttheme_minimal()
  panel <- ggarrange(plot_1 + theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "in")), 
                     plot_2 + theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "in")),
                     tableGrob(anova_df, theme = Tmin), 
                     tableGrob(random_effects_df, theme = Tmin, rows = NULL), 
                     tableGrob(table, theme = Tmin, rows = NULL), 
                     ncol = 1, heights = c(0.6, 0.6, 0.1, 0.2, 0.2),
                     labels = c("  Standard Boxplot ", "Predicted Mean ± 2*SE", "    Anova    ", "Random Effects", "   Contrasts  "))
  panel <- annotate_figure(panel, top = text_grob(
    paste0("Analysis for ", each_temp, "C timepoint ", each_endpoint, " h. P-value: ", cutoff_pvalue, " and FC: ", cutoff_FC),
                                                         face = "bold", size = 14, color = "red"))
  
  # Save files
  ggsave(panel, filename = paste0(figures_folder, "/summaryLDH_", each_temp, "C_", each_endpoint, "h.pdf"), width = 9, height = 14)
  ggsave(plot_2, filename = paste0(figures_folder, "/plotLDH_", each_temp, "C_", each_endpoint, "h.png"), width = 7, height = 6)
  saveRDS(plot_2, file.path(figures_folder, paste0("plotLDH_", each_temp, "C_", each_endpoint, "h.rds")))
  write_csv(anova_df, file.path(stats_folder, paste0("anova_", each_temp, "C_", each_endpoint, "h.csv")))
  write_csv(random_effects_df, file.path(stats_folder, paste0("stats_random_effects_", each_temp, "C_", each_endpoint, "h.csv")))
  write_csv(contrasts_df, file.path(stats_folder, paste0("stats_contrasts_", each_temp, "C_", each_endpoint, "h.csv")))
  write_csv(data_summary_df, file.path(stats_folder, paste0("stats_summary_", each_temp, "C_", each_endpoint, "h.csv")))
  
  return(list(
    anova = anova_df,
    random_effects = random_effects_df,
    contrasts = contrasts_df,
    data_summary = data_summary_df,
    data_subset = data_subset,
    plot_1 = plot_1,
    plot_2 = plot_2,
    panel = panel
  ))
}

Apply to each temp

Main Data: 34C

analysis_function(LDH_FC, each_temp = "34", each_endpoint = "48", cutoff_pvalue = 0.05, cutoff_FC = 1)
$anova
         Sum Sq Mean Sq NumDF   DenDF F value  Pr(>F) sign
Species 1.9e+02 6.3e+01     3 6.6e+01   4e+00 1.1e-02    *

$random_effects
        grp        var1 var2       vcov     sdcor proportion
1 Line:Date (Intercept) <NA>  0.0000000 0.0000000       0.00
2      Line (Intercept) <NA>  0.3787548 0.6154306       2.36
3  Residual        <NA> <NA> 15.6531696 3.9564087      97.64

$contrasts
           contrast condition1  condition2   estimate       SE       df
1 Dpi - Uncolonized        Dpi Uncolonized  0.5376209 1.261251 51.62153
2 Sau - Uncolonized        Sau Uncolonized  4.0243686 1.312822 53.38139
3 Spn - Uncolonized        Spn Uncolonized -0.2551308 1.289533 52.44458
     t.ratio     p.value p.adj.holm sign
1  0.4262599 0.671690735 1.00000000     
2  3.0654325 0.003406184 0.01021855    *
3 -0.1978474 0.843929327 1.00000000     
                                                 group1
1 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>
2 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>
3 <i><b><span style='color:#927ED1;'>Spn</span></i></b>
                                            group2 mean.predval.1
1 <b><span style='color:#5b5b5b;'>Uncol</span></b>       2.147787
2 <b><span style='color:#5b5b5b;'>Uncol</span></b>       5.634535
3 <b><span style='color:#5b5b5b;'>Uncol</span></b>       1.355036
  mean.predval.2        FC highlighted
1       1.610167  1.333891           +
2       1.610167  3.499349           +
3       1.610167 -1.188283           -

$data_summary
# A tibble: 4 × 7
# Groups:   Species [4]
  Species    bacteria_label mean.real mean.predval mean.predval.se   max     min
  <fct>      <fct>              <dbl>        <dbl>           <dbl> <dbl>   <dbl>
1 Dpi        <i><b><span s…      2.11         2.15           1.03   4.20  0.0973
2 Sau        <i><b><span s…      5.62         5.63           1.08   7.80  3.47  
3 Spn        <i><b><span s…      1.33         1.36           1.06   3.47 -0.758 
4 Uncoloniz… <b><span styl…      1.59         1.61           0.884  3.38 -0.158 

$data_subset
         Date   Line     Species
1  01_17_2023 HNO204         Dpi
2  01_17_2023 HNO204         Sau
3  01_17_2023 HNO204         Spn
4  01_17_2023 HNO204 Uncolonized
5  01_24_2023 HNO919         Dpi
6  01_24_2023 HNO919         Spn
7  01_24_2023 HNO919 Uncolonized
8  02_07_2023 HNO204         Dpi
9  02_07_2023 HNO204         Sau
10 02_07_2023 HNO204         Spn
11 02_07_2023 HNO204 Uncolonized
12 02_14_2023 HNO919         Dpi
13 02_14_2023 HNO919         Sau
14 02_14_2023 HNO919         Spn
15 02_14_2023 HNO919 Uncolonized
16 02_28_2023 HNO204         Dpi
17 02_28_2023 HNO204         Sau
18 02_28_2023 HNO204         Spn
19 02_28_2023 HNO204 Uncolonized
20 08_09_2022 HNO918         Dpi
21 08_09_2022 HNO918         Spn
22 08_09_2022 HNO918 Uncolonized
23 08_10_2022 HNO204         Dpi
24 08_10_2022 HNO204         Spn
25 08_10_2022 HNO204 Uncolonized
26 08_16_2022 HNO919         Dpi
27 08_16_2022 HNO919         Sau
28 08_16_2022 HNO919 Uncolonized
29 08_17_2022 HNO918         Dpi
30 08_17_2022 HNO918         Sau
31 08_17_2022 HNO918 Uncolonized
32 08_23_2022 HNO204         Dpi
33 08_23_2022 HNO204         Spn
34 08_23_2022 HNO204 Uncolonized
35 08_30_2022 HNO918         Dpi
36 08_30_2022 HNO918         Sau
37 08_30_2022 HNO918         Spn
38 08_30_2022 HNO918 Uncolonized
39 08_31_2022 HNO204         Dpi
40 08_31_2022 HNO204         Spn
41 08_31_2022 HNO204 Uncolonized
42 09_06_2022 HNO919         Dpi
43 09_06_2022 HNO919         Sau
44 09_06_2022 HNO919         Spn
45 09_06_2022 HNO919 Uncolonized
46 09_07_2022 HNO918         Dpi
47 09_07_2022 HNO918 Uncolonized
48 09_13_2022 HNO204         Dpi
49 09_13_2022 HNO204         Spn
50 09_13_2022 HNO204 Uncolonized
51 09_14_2022 HNO919         Dpi
52 09_14_2022 HNO919         Sau
53 09_14_2022 HNO919         Spn
54 09_14_2022 HNO919 Uncolonized
55 09_20_2022 HNO918         Sau
56 09_20_2022 HNO918 Uncolonized
57 09_21_2022 HNO204         Sau
58 09_21_2022 HNO204 Uncolonized
59 09_27_2023 HNO919         Sau
60 09_27_2023 HNO919         Spn
61 09_27_2023 HNO919 Uncolonized
62 10_04_2022 HNO918         Sau
63 10_04_2022 HNO918         Spn
64 10_04_2022 HNO918 Uncolonized
65 10_11_2022 HNO204         Sau
66 10_11_2022 HNO204         Spn
67 10_11_2022 HNO204 Uncolonized
68 10_25_2022 HNO204         Sau
69 10_25_2022 HNO204 Uncolonized
70 11_01_2022 HNO919         Dpi
71 11_01_2022 HNO919 Uncolonized
72 11_02_2022 HNO919 Uncolonized
                                          bacteria_label Temp Well_Endpoint
1  <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
2  <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
3  <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
4       <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
5  <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
6  <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
7       <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
8  <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
9  <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
10 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
11      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
12 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
13 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
14 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
15      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
16 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
17 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
18 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
19      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
20 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
21 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
22      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
23 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
24 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
25      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
26 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
27 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
28      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
29 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
30 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
31      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
32 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
33 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
34      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
35 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
36 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
37 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
38      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
39 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
40 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
41      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
42 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
43 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
44 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
45      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
46 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
47      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
48 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
49 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
50      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
51 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
52 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
53 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
54      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
55 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
56      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
57 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
58      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
59 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
60 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
61      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
62 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
63 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
64      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
65 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
66 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   34            48
67      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
68 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   34            48
69      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
70 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   34            48
71      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
72      <b><span style='color:#5b5b5b;'>Uncol</span></b>   34            48
           FC predval.fit predval.se.fit
1   1.6107056    2.147787      1.0252266
2   1.1521253    5.634535      1.0825826
3   1.2196382    1.355036      1.0566856
4   1.4673629    1.610167      0.8842437
5   1.1666667    2.147787      1.0252266
6   1.0420561    1.355036      1.0566856
7   1.0629921    1.610167      0.8842437
8   1.1233244    2.147787      1.0252266
9  18.7660099    5.634535      1.0825826
10  1.1213720    1.355036      1.0566856
11  1.3569554    1.610167      0.8842437
12  1.1079295    2.147787      1.0252266
13  1.4253394    5.634535      1.0825826
14  1.1241535    1.355036      1.0566856
15  0.8728324    1.610167      0.8842437
16  0.9670588    2.147787      1.0252266
17  1.1798715    5.634535      1.0825826
18  1.0730594    1.355036      1.0566856
19  1.2884097    1.610167      0.8842437
20  2.8243512    2.147787      1.0252266
21  1.7787115    1.355036      1.0566856
22  2.6411765    1.610167      0.8842437
23  2.3894340    2.147787      1.0252266
24  1.8271437    1.355036      1.0566856
25  2.0180878    1.610167      0.8842437
26  2.1775535    2.147787      1.0252266
27  1.7889578    5.634535      1.0825826
28  1.7968385    1.610167      0.8842437
29  2.4007134    2.147787      1.0252266
30  4.8096318    5.634535      1.0825826
31  1.8421673    1.610167      0.8842437
32  5.1042850    2.147787      1.0252266
33  0.9423572    1.355036      1.0566856
34  2.4878412    1.610167      0.8842437
35  2.8830189    2.147787      1.0252266
36  1.5291971    5.634535      1.0825826
37  2.1863799    1.355036      1.0566856
38  2.5090253    1.610167      0.8842437
39  3.6423077    2.147787      1.0252266
40  1.6666667    1.355036      1.0566856
41  1.6231884    1.610167      0.8842437
42  1.9311741    2.147787      1.0252266
43  1.9799197    5.634535      1.0825826
44  1.1392405    1.355036      1.0566856
45  1.2775330    1.610167      0.8842437
46  1.3724696    2.147787      1.0252266
47  1.7096774    1.610167      0.8842437
48  1.5791506    2.147787      1.0252266
49  1.3440000    1.355036      1.0566856
50  1.3512397    1.610167      0.8842437
51  2.3795380    2.147787      1.0252266
52  1.4391691    5.634535      1.0825826
53  1.1746988    1.355036      1.0566856
54  1.2951807    1.610167      0.8842437
55 30.6586103    5.634535      1.0825826
56  1.6335404    1.610167      0.8842437
57  1.2132353    5.634535      1.0825826
58  1.5528701    1.610167      0.8842437
59  1.8671679    5.634535      1.0825826
60  1.4087302    1.355036      1.0566856
61  1.6336634    1.610167      0.8842437
62  1.0116618    5.634535      1.0825826
63  0.9728097    1.355036      1.0566856
64  1.2942857    1.610167      0.8842437
65 13.7467018    5.634535      1.0825826
66  1.2047872    1.355036      1.0566856
67  1.3798701    1.610167      0.8842437
68  1.6891026    5.634535      1.0825826
69  1.6053333    1.610167      0.8842437
70  1.2621723    2.147787      1.0252266
71  1.3215686    1.610167      0.8842437
72  1.0997067    1.610167      0.8842437

$plot_1


$plot_2


$panel

Supplemental Data: 37C

analysis_function(LDH_FC, each_temp = "37", each_endpoint = "48", cutoff_pvalue = 0.05, cutoff_FC = 1)
$anova
         Sum Sq Mean Sq NumDF   DenDF F value  Pr(>F) sign
Species 3.4e+01 1.1e+01     3 3.2e+01 4.2e+00 1.3e-02    *

$random_effects
        grp        var1 var2      vcov     sdcor proportion
1 Line:Date (Intercept) <NA> 0.1592123 0.3990142       5.63
2      Line (Intercept) <NA> 0.0000000 0.0000000       0.00
3  Residual        <NA> <NA> 2.6703268 1.6341135      94.37

$contrasts
           contrast condition1  condition2   estimate        SE       df
1 Dpi - Uncolonized        Dpi Uncolonized  0.6737338 0.6154067 36.63814
2 Sau - Uncolonized        Sau Uncolonized  1.6196034 0.6319761 37.41009
3 Spn - Uncolonized        Spn Uncolonized -0.7311755 0.6330511 37.41164
    t.ratio    p.value p.adj.holm sign
1  1.094778 0.28075888 0.51082708     
2  2.562761 0.01453752 0.04361255    *
3 -1.155002 0.25541354 0.51082708     
                                                 group1
1 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>
2 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>
3 <i><b><span style='color:#927ED1;'>Spn</span></i></b>
                                            group2 mean.predval.1
1 <b><span style='color:#5b5b5b;'>Uncol</span></b>       2.681271
2 <b><span style='color:#5b5b5b;'>Uncol</span></b>       3.627141
3 <b><span style='color:#5b5b5b;'>Uncol</span></b>       1.276362
  mean.predval.2        FC highlighted
1       2.007537  1.335602           +
2       2.007537  1.806761           +
3       2.007537 -1.572859           -

$data_summary
# A tibble: 4 × 7
# Groups:   Species [4]
  Species     bacteria_label  mean.real mean.predval mean.predval.se   max   min
  <fct>       <fct>               <dbl>        <dbl>           <dbl> <dbl> <dbl>
1 Dpi         <i><b><span st…      2.68         2.68           0.485  3.65 1.71 
2 Sau         <i><b><span st…      3.60         3.63           0.507  4.64 2.61 
3 Spn         <i><b><span st…      1.28         1.28           0.506  2.29 0.263
4 Uncolonized <b><span style…      2.01         2.01           0.386  2.78 1.24 

$data_subset
         Date   Line     Species
1  02_28_2023 HNO204         Sau
2  02_28_2023 HNO204 Uncolonized
3  08_09_2022 HNO918         Dpi
4  08_09_2022 HNO918         Spn
5  08_09_2022 HNO918 Uncolonized
6  08_10_2022 HNO204         Dpi
7  08_10_2022 HNO204         Spn
8  08_10_2022 HNO204 Uncolonized
9  08_17_2022 HNO918         Dpi
10 08_17_2022 HNO918         Sau
11 08_17_2022 HNO918         Spn
12 08_17_2022 HNO918 Uncolonized
13 08_23_2022 HNO204         Dpi
14 08_23_2022 HNO204         Spn
15 08_23_2022 HNO204 Uncolonized
16 08_30_2022 HNO918         Dpi
17 08_30_2022 HNO918         Sau
18 08_30_2022 HNO918         Spn
19 08_30_2022 HNO918 Uncolonized
20 08_31_2022 HNO204         Dpi
21 08_31_2022 HNO204         Spn
22 08_31_2022 HNO204 Uncolonized
23 09_06_2022 HNO919         Dpi
24 09_06_2022 HNO919         Sau
25 09_06_2022 HNO919         Spn
26 09_06_2022 HNO919 Uncolonized
27 09_07_2022 HNO918         Dpi
28 09_07_2022 HNO918 Uncolonized
29 09_13_2022 HNO204         Dpi
30 09_13_2022 HNO204         Spn
31 09_13_2022 HNO204 Uncolonized
32 09_14_2022 HNO919         Dpi
33 09_14_2022 HNO919         Sau
34 09_14_2022 HNO919 Uncolonized
35 09_20_2022 HNO918         Sau
36 09_20_2022 HNO918 Uncolonized
37 09_21_2022 HNO204         Sau
38 09_21_2022 HNO204 Uncolonized
39 09_27_2023 HNO919         Sau
40 09_27_2023 HNO919         Spn
41 09_27_2023 HNO919 Uncolonized
42 10_04_2022 HNO918         Sau
43 10_04_2022 HNO918         Spn
44 10_04_2022 HNO918 Uncolonized
45 10_11_2022 HNO204         Sau
46 10_11_2022 HNO204         Spn
47 10_11_2022 HNO204 Uncolonized
48 10_25_2022 HNO204         Sau
49 10_25_2022 HNO204 Uncolonized
50 11_01_2022 HNO919         Dpi
51 11_01_2022 HNO919 Uncolonized
52 11_02_2022 HNO919         Dpi
53 11_02_2022 HNO919 Uncolonized
                                          bacteria_label Temp Well_Endpoint
1  <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
2       <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
3  <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
4  <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
5       <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
6  <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
7  <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
8       <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
9  <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
10 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
11 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
12      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
13 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
14 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
15      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
16 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
17 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
18 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
19      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
20 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
21 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
22      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
23 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
24 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
25 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
26      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
27 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
28      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
29 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
30 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
31      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
32 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
33 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
34      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
35 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
36      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
37 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
38      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
39 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
40 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
41      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
42 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
43 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
44      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
45 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
46 <i><b><span style='color:#927ED1;'>Spn</span></i></b>   37            48
47      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
48 <i><b><span style='color:#6D05A0;'>Sau</span></i></b>   37            48
49      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
50 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
51      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
52 <i><b><span style='color:#0443DC;'>Dpi</span></i></b>   37            48
53      <b><span style='color:#5b5b5b;'>Uncol</span></b>   37            48
           FC predval.fit predval.se.fit
1   1.5136986    3.627141      0.5065006
2   1.4781705    2.007537      0.3859056
3   2.2846571    2.681271      0.4850404
4   1.4300518    1.276362      0.5064943
5   2.5632816    2.007537      0.3859056
6   2.7198708    2.681271      0.4850404
7   0.2867798    1.276362      0.5064943
8   2.5420207    2.007537      0.3859056
9   3.0947136    2.681271      0.4850404
10  2.2683449    3.627141      0.5065006
11  0.4345635    1.276362      0.5064943
12  1.7299169    2.007537      0.3859056
13  2.2388127    2.681271      0.4850404
14  1.7504621    1.276362      0.5064943
15  3.1889793    2.007537      0.3859056
16  2.9825175    2.681271      0.4850404
17  3.1650165    3.627141      0.5065006
18  1.9658703    1.276362      0.5064943
19  2.1770833    2.007537      0.3859056
20  3.2170819    2.681271      0.4850404
21  1.1304348    1.276362      0.5064943
22  4.0029940    2.007537      0.3859056
23  2.4717742    2.681271      0.4850404
24  1.8406375    3.627141      0.5065006
25  1.2559055    1.276362      0.5064943
26  1.3755102    2.007537      0.3859056
27  2.7200000    2.681271      0.4850404
28  3.1715116    2.007537      0.3859056
29  4.2857143    2.681271      0.4850404
30  1.6109589    1.276362      0.5064943
31  2.2633745    2.007537      0.3859056
32  2.5892473    2.681271      0.4850404
33  3.8017429    3.627141      0.5065006
34  0.8404352    2.007537      0.3859056
35  1.5977444    3.627141      0.5065006
36  1.5537190    2.007537      0.3859056
37  1.6948640    3.627141      0.5065006
38  1.2089041    2.007537      0.3859056
39  1.8498728    3.627141      0.5065006
40  1.3581081    1.276362      0.5064943
41  0.9537037    2.007537      0.3859056
42  7.2729885    3.627141      0.5065006
43  1.2089552    1.276362      0.5064943
44  1.4824561    2.007537      0.3859056
45  1.9615385    3.627141      0.5065006
46  1.6141975    1.276362      0.5064943
47  1.2881844    2.007537      0.3859056
48 12.6428571    3.627141      0.5065006
49  1.7321937    2.007537      0.3859056
50  1.7596899    2.681271      0.4850404
51  2.4734848    2.007537      0.3859056
52  1.8522727    2.681271      0.4850404
53  2.1172840    2.007537      0.3859056

$plot_1


$plot_2


$panel

Merged Files

pdf_output <- file.path(figures_folder, "HNOBac_SummaryLDH.pdf")

# Check if the file exists, and delete if it does
if (file.exists(pdf_output)) {
  file.remove(pdf_output)
}
[1] TRUE
# Now combine the PDFs into a new file
pdf_files <- list.files(figures_folder, pattern = "\\.pdf$", full.names = TRUE)
qpdf::pdf_combine(input = pdf_files, output = pdf_output)
[1] "data/outputs/LDH/figures/HNOBac_SummaryLDH.pdf"
CFUs
Cytokines
Source Code
---
execute:
  message: FALSE
  warning: FALSE
---

# LDH {.unnumbered}

```{r}
library(tidyverse)
library(readxl)
library(ggpubr)
library(ggtext)
library(scales)
library(lme4)
library(afex)
library(emmeans)
library(gridExtra)
```

## Data Input and Selection

### File Paths

```{r}
# Folder paths
input_path <- "data/input_data/LDH/"
metadata_path <- "data/metadata/LDH"

# Create subfolders for output files
dataframes_folder <- "data/dataframes"
if (!file.exists("data/dataframes")) {
  dir.create("data/dataframes", recursive = TRUE)
}

# Load data and metadata
input_data <- read_excel(file.path(input_path, "LAK24_0118_LDHData.xlsx")) %>% mutate_if(is.character, factor)
Bac_order <- read_csv(file.path(metadata_path, "Order_Bacteria_LDH.csv"))
```

### Data clean-up

```{r}
# Filter only samples matched to CFU wells
LDH_data <- input_data %>% 
  filter(CFUmatched == "T")

# Factor Ordering and Styling
LDH_data <- merge(LDH_data, Bac_order, by = "Species") 
LDH_data$bacteria_label <- factor(LDH_data$bacteria_label, levels = Bac_order$bacteria_label)
LDH_data$Line <- fct_recode(LDH_data$Line, "HNO918" = "A", "HNO204" = "B", "HNO919" = "C") 

# Average technical replicates for each individual experiment
LDH_avg <- LDH_data %>%
  group_by(Date, Line, Species, bacteria_label, Temp, Well_Endpoint, Collection_Time) %>% 
  summarise(avg_Value = mean(Value)) %>%
  # Add labels to distinguish final collection times from time 0
  mutate(Collection_Label = if_else(Collection_Time == 24 | Collection_Time == 48, "final", "initial"))

# Calculate log2 fold change between final time point vs. time 0
LDH_FC <- LDH_avg %>%
  group_by(Date, Line, Species, bacteria_label, Temp, Well_Endpoint) %>%
  reframe(FC = avg_Value[Collection_Label == "final"]/avg_Value[Collection_Label == "initial"]) 
```

### Saving files

```{r}
# Save data frames as CSV files in the output folder
write_csv(LDH_data, file.path(dataframes_folder, "LDH_values.csv"))
write_csv(LDH_FC, file.path(dataframes_folder, "LDH_FC.csv"))

# Save data frames as R objects in the output folder
saveRDS(LDH_data, file.path(dataframes_folder, "LDH_values.rds"))
saveRDS(LDH_FC, file.path(dataframes_folder, "LDH_FC.rds"))

# Cleaning-up all objects from the environment
rm(list = ls())

# Use this to read the final objects
LDH_data <- readRDS("data/dataframes/LDH_values.rds")
LDH_FC <- readRDS("data/dataframes/LDH_FC.rds")
```

## Stats and Plots

### File Paths

```{r}
# Folder paths
dataframes_path <- "data/dataframes"
metadata_path <- "data/metadata/LDH"

# Create subfolders for output files
figures_folder <- "data/outputs/LDH/figures"
if (!file.exists("data/outputs/LDH/figures")) {
  dir.create("data/outputs/LDH/figures", recursive = TRUE)
}
stats_folder <- "data/outputs/LDH/stats"
if (!file.exists("data/outputs/LDH/stats")) {
  dir.create("data/outputs/LDH/stats", recursive = TRUE)
}

# Load data and metadata
LDH_FC <- readRDS("data/dataframes/LDH_FC.rds")
Bac_order <- read_csv(file.path(metadata_path, "Order_Bacteria_LDH.csv"))
```

### Function for each temp condition and endpoint

```{r}
# Function to analyze each temp condition
analysis_function <- function(data, each_temp, each_endpoint, cutoff_pvalue, cutoff_FC) {
  
  # Subset the data to the selected temp
  data_subset <- data %>%
    filter(Temp == each_temp) %>%
    filter(Well_Endpoint == each_endpoint)
  
  # Mixed-effects model with random effects
  model <- lmer(FC ~ Species 
                + (1|Line) + (1|Line:Date), 
                data = data_subset)
  #Anova
  anova <- anova(model)
  anova_df <- as.data.frame(anova) %>%
    mutate(sign = case_when(
      `Pr(>F)` < cutoff_pvalue ~ "*",
      TRUE ~ "")) %>%
    mutate_if(is.numeric, ~ format(., digits = 2, scientific = TRUE))
  
  # Calculate Individual contrasts
  emmeans_model <- emmeans(model, ~ Species)
  emmeans_contrasts <- pairs(emmeans_model, adjust = "none")    
  
  # Extract random effects and convert to dataframe (if not singular)
  random_effects_df <- as.data.frame(VarCorr(model)) %>%
    mutate(proportion = round(100 * (vcov / sum(vcov)), 2)) 
  
  # Adds predictions based on fixed effects, averaged over random effects. It gives a population estimate
  data_subset <- cbind(data_subset, predval = predict(model,re.form = NA, se.fit = TRUE))
  data_summary_df <- data_subset %>%
    group_by(Species, bacteria_label) %>%
    summarize(mean.real = mean(FC),
              mean.predval = mean(predval.fit), 
              mean.predval.se = mean(predval.se.fit)) %>%
    mutate(max = mean.predval + 2*mean.predval.se,
           min = mean.predval - 2*mean.predval.se)
  
  # Convert contrasts to dataframe and adjust pvalues. Filter contrast to Uncolonized only
  contrasts_df <- as.data.frame(summary(emmeans_contrasts)) %>%
    filter(str_detect(contrast, "Uncolonized")) %>%
    mutate(p.adj.holm = p.adjust(p.value, method = "holm")) %>%
    mutate(sign = case_when(
      p.adj.holm < 0.05 ~ "*",
      TRUE ~ ""))
  
  # Edits to the contrast dataframe to include pvalue brackets in plot
  contrasts_df <- contrasts_df %>%
    separate(contrast, into = c("condition1", "condition2"), sep = " - ", remove = F) 
  
  contrasts_df$group1 <- Bac_order$bacteria_label[match(contrasts_df$condition1, Bac_order$Species)]
  contrasts_df$group2 <- Bac_order$bacteria_label[match(contrasts_df$condition2, Bac_order$Species)]
  
  # Calculate fold-change values for each contrast
  contrasts_df <- contrasts_df %>%
    ungroup() %>%
    left_join(select(data_summary_df, Species, mean.predval), by = join_by(condition1 == Species)) %>%
    left_join(select(data_summary_df, Species, mean.predval), by = join_by(condition2 == Species), suffix = c(".1", ".2")) %>%
    mutate(FC = mean.predval.1 / mean.predval.2,
           FC = if_else(FC < 1, -1 / FC, FC),
           highlighted = case_when(
             FC <= -cutoff_FC ~ "-",
             FC >= cutoff_FC ~ "+",
             TRUE ~ "")) 
  
  # Select p values to plot and define their location
  contrast_sign <- contrasts_df %>%
    filter(sign != "" & highlighted != "") %>%
    mutate(p.adj.holm = format(p.adj.holm, digits = 2, scientific = TRUE))
  
  location <- max(data_subset$FC, na.rm = TRUE) * 1.1
  
  # Standard Boxplot
  plot_1 <- ggplot() +
    geom_boxplot(data = data_subset, aes(x = bacteria_label, y = FC, fill = bacteria_label)) + 
    
    geom_jitter(data = data_subset, aes(x = bacteria_label, y = FC, shape = Line), 
                fill = "grey50", color = "grey30", size = 2, width = 0.05, stroke = 0.75) +
    
    scale_fill_manual(values = c("grey80","#AA35E3","#2e67f2","#927ed1")) +
    scale_shape_manual(values = c(21, 22, 24)) +
    
    labs(x = "",
         y = "Fold change in LDH from -1 hour") +
    
    theme_bw() +
    theme(panel.grid = element_blank(),
          legend.position = "none",
          text = element_text(size = 20), 
          axis.text.x = element_markdown(), 
          axis.text.y = element_text(color = "black"))
  
  # Plot with predicted means and standard errors of the estimates
  plot_2 <- ggplot() +
    geom_point(data = data_subset, 
                aes(x = bacteria_label, y = FC, fill = bacteria_label, color = bacteria_label, shape = Line), 
                position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
                size = 1.5, alpha = 0.75, stroke = 0.75) +
    
    geom_point(data = data_summary_df, aes(x = bacteria_label, y = mean.predval), shape = 3, size = 3) +
    geom_errorbar(data = data_summary_df, aes(x = bacteria_label,
                                              y = mean.predval,
                                              ymax = max,
                                              ymin = min),
                  width = 0.4) +
    
    scale_fill_manual(values = c("#5b5b5b","#AA35E3","#2e67f2","#927ed1")) +
    scale_color_manual(values = c("#5b5b5b","#AA35E3","#2e67f2","#927ed1")) +
    scale_shape_manual(values = c(21, 22, 24)) +
    
    scale_y_continuous(expand = c(0.1,0)) +
    
    labs(x = "",
         y = "Fold change in LDH from -1 hour",
         fill = "Bacteria", color = "Bacteria", shape = "HNO Line") +
    
    theme_bw() +
    theme(panel.grid = element_blank(), 
          legend.text = element_markdown(),
          text = element_text(size = 20), 
          axis.text.x = element_markdown(), 
          axis.text.y = element_text(color = "black"))
  
  # Conditionally add p-value annotations layer
  if (nrow(contrast_sign) > 0) {
    plot_2 <- plot_2 +
      stat_pvalue_manual(contrast_sign, label = "p.adj.holm", y.position = location,
                         tip.length = 0.02, bracket.shorten = 0.2, vjust = -0.2, bracket.size = 0.3, size = 3.5)
  } else {
    plot_2 <- plot_2
  }
  
  # Arrange plot and tables for summary pdf
  table <- contrasts_df %>%
    select(condition1, condition2, p.adj.holm, sign, mean.predval.1, mean.predval.2, FC, highlighted) %>%
    mutate(p.adj.holm = format(p.adj.holm, digits = 2, scientific = TRUE))
  
  Tmin <- ttheme_minimal()
  panel <- ggarrange(plot_1 + theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "in")), 
                     plot_2 + theme(plot.margin = unit(c(0.25,0.25,0.25,0.25), "in")),
                     tableGrob(anova_df, theme = Tmin), 
                     tableGrob(random_effects_df, theme = Tmin, rows = NULL), 
                     tableGrob(table, theme = Tmin, rows = NULL), 
                     ncol = 1, heights = c(0.6, 0.6, 0.1, 0.2, 0.2),
                     labels = c("  Standard Boxplot ", "Predicted Mean ± 2*SE", "    Anova    ", "Random Effects", "   Contrasts  "))
  panel <- annotate_figure(panel, top = text_grob(
    paste0("Analysis for ", each_temp, "C timepoint ", each_endpoint, " h. P-value: ", cutoff_pvalue, " and FC: ", cutoff_FC),
                                                         face = "bold", size = 14, color = "red"))
  
  # Save files
  ggsave(panel, filename = paste0(figures_folder, "/summaryLDH_", each_temp, "C_", each_endpoint, "h.pdf"), width = 9, height = 14)
  ggsave(plot_2, filename = paste0(figures_folder, "/plotLDH_", each_temp, "C_", each_endpoint, "h.png"), width = 7, height = 6)
  saveRDS(plot_2, file.path(figures_folder, paste0("plotLDH_", each_temp, "C_", each_endpoint, "h.rds")))
  write_csv(anova_df, file.path(stats_folder, paste0("anova_", each_temp, "C_", each_endpoint, "h.csv")))
  write_csv(random_effects_df, file.path(stats_folder, paste0("stats_random_effects_", each_temp, "C_", each_endpoint, "h.csv")))
  write_csv(contrasts_df, file.path(stats_folder, paste0("stats_contrasts_", each_temp, "C_", each_endpoint, "h.csv")))
  write_csv(data_summary_df, file.path(stats_folder, paste0("stats_summary_", each_temp, "C_", each_endpoint, "h.csv")))
  
  return(list(
    anova = anova_df,
    random_effects = random_effects_df,
    contrasts = contrasts_df,
    data_summary = data_summary_df,
    data_subset = data_subset,
    plot_1 = plot_1,
    plot_2 = plot_2,
    panel = panel
  ))
}
```

### Apply to each temp

#### Main Data: 34C

```{r}
#| warning: FALSE
analysis_function(LDH_FC, each_temp = "34", each_endpoint = "48", cutoff_pvalue = 0.05, cutoff_FC = 1)
```

#### Supplemental Data: 37C

```{r}
#| warning: FALSE
analysis_function(LDH_FC, each_temp = "37", each_endpoint = "48", cutoff_pvalue = 0.05, cutoff_FC = 1)
```

#### Merged Files

```{r}
pdf_output <- file.path(figures_folder, "HNOBac_SummaryLDH.pdf")

# Check if the file exists, and delete if it does
if (file.exists(pdf_output)) {
  file.remove(pdf_output)
}

# Now combine the PDFs into a new file
pdf_files <- list.files(figures_folder, pattern = "\\.pdf$", full.names = TRUE)
qpdf::pdf_combine(input = pdf_files, output = pdf_output)
```